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ABSTRACT 

We describe a new finite element method (FEM) to construct continuous equilib- 
rium distribution functions of stellar Systems. The method is a generalization of 
Schwarzschild's orbit superposition method from the Space of discrete functions to 
continuous ones. In contrast to Schwarzschild's method, FEM produces a continuous 
distribution function (DF) and satisfies the intra element continuity and Jeans equa- 
tions. The method employs two finite-element meshes, one in configuration space and 
one in action space. The DF is represented by its values at the nodes of the action-space 
mesh and by interpolating functions inside the Clements. The Galerkin projection of 
all equations that involve the DF leads to a linear System of equations, which can be 
solved for the nodal values of the DF using linear or quadratic programming, or other 
optimization methods. We illustrate the superior Performance of FEM by construct- 
ing ergodic and anisotropic equilibrium DFs for spherical stellar Systems (Hcrnquist 
modeis). We also show that explicitly constraining the DF by the Jeans equations 
leads to smoother and/or more accurate Solutions with both Schwarzschild's method 
and FEM. 

Key words: stellar dynamics, methods: numerical, galaxies: kinematics and dynam- 
ics, galaxies: elliptical 



1 INTRODUCTION 

For over three decades, Schwarzschild's (1979) orbit super- 
position method has been one of the most important numer- 
ical tools for modelling the equilibrium stat es of spherical 
IIRichstone &: Tremainelll984l ). axisymmetric (|Thonias et al.l 
[2OOJ) and triaxial galaxies (van den Bosch et al. 2008 and 
references therein). Schwarzschild's method constructs dis- 
crete phase-space distribution functions (DF) and works 
even if t he gravitational potential S upports chaotic Or- 
bits (e.g.. lCapuzzo-Dolcetta et al.ll2007l V Observational con- 
straints can also be incorporated, in particular on the sur- 
face brightness or the line-of-sight velocity distributions. 
Schwarzschild's basic assumptions were: (i) the amount of 
mass contributed by an orbit to a cell/element in the con- 
figuration Space is proportional to the fraction of time spent 
by that orbit inside the element; (ii) the matter density and 
velocity distribution inside each element is constant; (iü) the 
DF is non-zero only on a subset of phase space with mea- 
sure zero (except perhaps in some cases where the potential 
admits large-scale chaos). In particular if all orbits are reg- 
ulär the DF is discrete, i.e., non-zero only at a finite set of 
positions in action space. 
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Despite its central role in modelling galaxies, Schwarz- 
schild's method has several shortcomings. (i) There is no 
mathematical proof that increasing the number of Clements 
and orbits in this scheine will guarantee the convergence of 
the coarse-grained DF to a smooth function. (ii) In practice, 
Schwarzschild modeis often converge rather slowly, in part 
because of the inverse square-root singularity in the density 
contributed by an orbit near a turning point. (iü) Working 
with discrete DFs is not ideal if we require their derivatives 
for linear stability analysis, or use them to set initial condi- 
tions for A^-body simulations. 

The majority of these limitations can be removed if one 
extends the method from the space of discrete functions to a 
more general continuous class. This is a problem in galactic 
dynamics whose Solution is overdue and we aim to solve it 
using a modified finite element method (FEM). 

The mathematical principles of FEM are rather simple. 
Assume a general governing equation 



A[u{x,t)] =0, 



(1) 



for the physical quantity u{x, t) with A being a partial 
integro-differential Operator, and seek the Solutions in terms 
of the coordinates x and the time t. An approximate Solution 
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of (Hl may be expanded as the series 

Jmax 

«(a;,t)= ^i*,(i)^,-(x), (2) 

where tpj{x) constitute a complete set of basis functions that 
satisfy any given boundary conditions. Substituting from ((2|) 
into m, multiplying both sides of the resulting equation by 
tpji[x), and integrating over the x doniain resuhs in 

/ 1pj, {x) A[u{x, t)\ AX =0, j' = 1, 2, . . . , jniax- (3) 

This is called the weighted residual form, weak form, or 
Galerkin projection of the governing equation U]), from 
which one may compute the unknown functions Uj{t). Here, 
the basic assumption is that if the equations (|3]) are satisfied 
for aU j' then equation ([T]) is also satisfied to an adequate 
degree of approximation. 

It is, however, a non-trivial and sometimes impossible 
task to find suitable basis functions that (i) satisfy all bound- 
ary conditions, (ii) make a complete set, (üi) do not con- 
tribute noise when the Solution is changing rapidly. More- 
over, the Integration over the spatial variables is expensive if 
the domain of each basis function is the entire x-space. One 
can overcome these difüculties by dividing the a;-space into 
A'^ finite elements. The union of elements, each of volume 
Vn, is the entire x domain. Each element contains A'^d nodes 
on its boundaries or in its inferior. The function u is approx- 
imated inside the nth element as a weighted sum of smooth 
shape functions ipnj{x); these are defined only within ele- 
ment n and are zero at all nodes except node j of bin n, 
where they are unity. The weights Unj are time-dependent 
and are chosen to fit the unknown value of u at the nodes. 
The determining equations of the weights are 



i>„j{x) A[u{x,t)] Ax 



ij'njix) A 



^U„j{t)lp„j{x) 



dx=0. 



(4) 



n = l,2,...,iV, j 



,Ni 



This procedure, which leaves us with a System of 
ordinary integro-differential equations for the weights 
Ujkjt) (the nodal values of u) is called the FEM 
ijZienkiewicz. Tavlor fc Zhull2005l l. 

For example, when the Operator A has the form Au = 
du/dt+Cu with jC{t) a linear Operator, the weighted residual 
form takes a simple matricial form 






A(i) • u(i) + F{t), 



(5) 



where the matrix A(t) is the projection of the Opera- 
tor —C and F{t) is a forcing vector. The vector u{t) 
contains nodal values of u{x,t). Therefore, the combina- 
tion of finite elements and Galerkin projection reduces an 
infinite-dimensional partial differential equation to a finite- 
dimensional System. 

A formula tion of FEM for stellar Systems was presented 
in ljalalil l|20ld ). where the perturbed coUisionless Boltzmann 
equation (CBE) was reduced to a form like ((5]) and solved 
over a ränge of finite ring elements in the configuration 
Space. That analysis, however, cannot be directly used to 
construct equilibrium DFs because they are not unique: ac- 



cording to the Jeans theorem, any distribution function / 
that depends on phase-space coordinates only through the 
Integrals of motion /is an equilibrium Solution of the CBE. 
Consequently the local Variation of / in the /-space is free. 
This implies non-uniqueness and / can admit discrete, piece- 
wise continuous, continuous, and differentiable Solutions. 

In this paper we develop a general method to build nu- 
merical DFs that exhibit nice properties of local differen- 
tiability and global continuity. After defining the problem 
in i jl.ll in !j2]we discuss finite elements, Interpolation func- 
tions, and their properties both in the configuration and 
action Spaces. In iJ31 we obtain the Galerkin projections of 
velocity moments, and the continuity and Jeans equations. 
Schwarzschild's method is derived as a special case of FEM 
in iJ3] and finite element modeis of spherical Systems are dis- 
cussed in [JS] We apply the FEM to the spherical Hernquist 
model in 36] and 37] contains a discussion of our results. 



1.1 Equilibrium stellar Systems 

The DF of a coUisionless System in dynamical equilibrium 
depends on the position x = {x-i_,X2,X3) and the velocity 
V = (vi, V2,V3) vectors only through the Integrals of motion. 

Integrable Systems The Hamilton-Jacobi equation is sep- 
arable in spherical Systems, razor-thin axisymmetric discs, 
and triaxial Systems in which the potential is of Stäckel form. 
Orbits in these Systems are regulär, and can be represented 
using a suitable action vector J — (Ji, J2, J3) and its asso- 
ciated angle variables w — {'Wi,W2,W3). The Hamiltonian 
function "H depends only on the action vector and the evo- 
lution of the angle variables is linear in time: 

w{t)^nt + w{o), n{j) = -^, (6) 

where fi = (^1,^2,^73) is the vector of orbital frequencies. 
The actions are isolating Integrals of motion, so any DF 
of the form /( J) defines a possible equilibrium stellar Sys- 
tem. For separable Systems, the actions can be computed by 
quadratures. In this study, we focus on building equilibrium 
modeis of stellar Systems with integrable potentials. 

Non-integrable Systems Generic axisymmetric or triax- 
ial Potentials are not integrable. There are surviving invari- 
ant fori of regulär orbits (cf. KAM theory) but these are sep- 
arated by chaotic layers. If the chaotic layers are narrow (the 
Potential is 'nearly integrable') the DF may be assumed to 
be zero in the chaotic phase subspace and written as a func- 
tion of the actions in the regulär phase subspace. However, 
the actions must then be calculated either using canonical 
perturbatio n theory or by generating a frequency map of 
the S ystem (JLaskarJllQQÖI : lnunterlbood iBinnev fc Tremain^ 
I2OO8I ). 



A simpler and more powerful approach, which can be 
used even if the potential is far from integrability, is to ex- 
press t he Integrals of motion in terms of initial conditions of 
orbits (jSchwarzschildl 1 1979 ) . Thus let [xo{x, v) , vo{x, v)] be 
the Position and velocity of the trajectory through (x, v) on 
some global surface of section T> through which all orbits 
must pass (e.g., a symmetry plane of a triaxial potential). 
Then any DF of the form f{xo,vo) defines an equilibrium 
stellar System. 
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1.1.1 Moments of the DF 

A coUisionless System with a given density function p{x) is 
a possible equilibrium if one can find a DF /( J) > so that 



Pix) = / /( J) dv. 



(7) 



We shall also sometimes use the first- and second-order ve- 
locity moments; 



u{x) = p{vi) (x) = Vi f{J) dv, 
t''^{x) = p{viVj){x)= / ViVjf{J) dv. 



(8) 
(9) 



In an equilibrium System these are related by the steady- 
state continuity equation 



^ OXi 

z — 1 

and Jeans equations 

3 

dxi '^ dx 



l^l^ = -Pn;r' ^ = 1,2,3, 



(10) 



(11) 



j=i 



with $(a;) being the potential. In Systems with spherical 
symmetry only the radial and tangential velocity dispersions 
matter and three Jeans equations reduce to one. 

We shall argue below that including constraints based 
on the continuity and Jeans equations can significantly im- 
prove the accuracy of both Schwarzschild and FEM modeis 
of stellar Systems. 



2 FINITE ELEMENTS IN CONFIGURATION 
AND ACTION SPACE 

We assume that the configuration space has been split into 
A'^ elements, each of A'^d nodes, and that the density p{x) 
is known at the nodal points. Inside each element, the den- 
sity function can be approximated by suitable interpolation 
(shape) functions. Denoting p„(x) as the functional form of 
the density inside the nth element, one may write 



p{x) = '^Hn{x) p„{x), pn{x) = '^gk,n{x) Pk,: 



(12) 



The density at the fcth node of the nth element has been 
indexed by the pair (fc,n). The function Hn{x) is unity in- 
side the nth element in the a;-space and zero outside. The 
interpolation functions gj^n{x) have the following properties: 



gj.niXkn) =5jk, j,k =1,2,. 



,Na 



(13) 



where Sjk is the Kronecker delta, and Xkn is the position vec- 
tor of the fcth node of the nth element. Figure[T]shows some 
elementary one-, two- and three-dimensional elements. The 
rectangular and brick elements can be distorted to obtain 
the s o-called mapped elements (jZienkiewicz. Taylor fc Zhul 
I2OO5I ). which help to reconstruct complex morphologies in 
curvilinear coordinates. For instance, elements confined be- 
tween confocal ellipsoids and hyperboloids can better de- 
scribe elliptical galaxy modeis that may have Potentials close 
to Stäckel form. Thin rings and spherical shells are the most 




2D rectangular 



^ 



V 



3D brick 



Figure 1. Elementary finite elements. The degree of interpolating 
polynomial increases by adding interior nodes, which can lie on 
the edges, sides or even inside elements. 



efficient elements for axisymmetric discs and spherical Sys- 
tems, respectively. 

Using the superscript T to transpose a vector/matrix, 
we define the row vector 

9n{x) = [ ffi,n(a;) ff2,n(a;) ••• 5iiv^,„(a;) ] , (14) 
and the column vector 



[ Pl,n p2,' 



PJVd 



(15) 



and rewrite the components of p„{x) in the following com- 
pact form 



Pn = 3n ■ bu. 



(16) 



Here a dot denotes matrix/vector multiplication. The above 
procedure can be readily applied to higher order velocity 
moments. In particular, we obtain 



i i ij jij 



(17) 



where the column vectors cjj and d'^ contain, respectively, 
the nodal values of u^ and r'-' inside the nth element. 

To construct a DF /(J), we divide the action space to 
M finite elements, each of Md nodes, and write 



/(J) = Y. Hm{J) f^{J) ■ P„ 



(18) 



with f^{J) and p^ being Md dimensional row and column 
vectors, respectively. The elements of the interpolating row 
vector f^{J) are denoted by fj.m{J) and they satisfy the 
condition 



Jj,m \Jk' 



Ojk, 



j,k = 1,2, 



,Md 



(19) 



for the action vector Jkm associated with the feth node of 
the mth element in action space. The union of the domains 
of fj^rn Covers the action space. 

In this study we use interpolation functions of Co class 
both in the x and J Spaces. The use of Co functions im- 
plies that all physical quantities are smooth (continuous 
and differentiable) inside each element and along its bound- 
ary lines. In the direction perpendicular to the boundary 
lines and at the nodes, the DF, density and higher order ve- 
locity moments will only be continuous. For example, con- 
sider the simplest one-dimensional set of Co elements: ele- 
ment n has boundaries at x„ and Xn+i > Xn and has two 
nodes, with node 1 at the smaller boundary Xn and node 2 
at the larger. The continuity of p{x) ai x = x„+i implies 
Pn{xn+i) = pn+i{xn+i). Usiug (fTS)) . this coudition reduces 
to 



&2,7i = 61, („+1). 



(20) 
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The first derivative of pn{^) with respect to x exists inside 
element n and is given by 



dpn , dgi,n , , dg2,n 

-^ — ~ Öl,n— ;t h 02.n—^ , 

OX OX OX 



(21) 



but the differentiability condition at the nodes of elements 
is not necessarily satisfied, i.e., 



dp„ 



dx 



+ 



dp 



n + l 



dx 



(22) 



One can resolve this problem by applying Ci finite elements. 
The application of C\ elements requires larger vectors of 
nodal quantities (which should now include partial deriva- 
tives), and thus larger element matrices. In this paper we 
restrict ourselves to Co elements; however, we note that C\ 
elements provide smoother Solutions (at the cost of larger 
matrices and greater analytic complexity), and are hkely to 
be useful when the partial derivatives of /( J) are also present 
(e.g., in linear stability analyses). 

Any vectorial function of the form f / v^ H-g (x )g„ (x) 
can b e expressed in terms of angle-action variables IjJalalil 
I2OIOI ): 



'j'^v'/H„{x)g^{x) = ^gj,(i,j,Zi,/2,n, J)e'* 



(23) 



with gl = g_j,. Here the asterisk Stands for complex conju- 
gation, fc is a 3-vector of integers and i = \/— T. To simplify 
the notation, we will denote g^i^j, h, h^n, J) by gs,(n, J) if 
h = h = 0, by 'giJyi, n, J) \i I2 = Q and h = 1, and by 
g^{i,j,n,X) if h = h = 1. These special cases correspond 
to the zeroth-, first-, and second-order velocity moments. 
The row vector g^ has the same dimension as g„ and it is 
calculated from 



9k 



1 



(27r)^ 



'1 h 



H„{x)g^{x) e ' "" dm 



(24) 



When a fest particle with the action vector J is inside the 
nth element in the configuration space, the function H„{x) is 
unity and that particle contributes to g,,. In other situations, 
the integrand of (I24|l will vanish. 

To compute g^,, we simply integrale the equations of 
motion corresponding to the action J or the initial con- 
ditions (a^j, «0) until the particle enters the nth element 
at time ti,„ and exits at t2,n- We then calculate the val- 
ues of the angles at the entry and exit times, wi^n and 
W2,n = wi,n + n(i2,n — ti,n)- Wc theu Carry out the Inte- 
gration (I24|l using Gaussian quadrature, typically with 8-15 
points. The numerical Integration of the equations of motion 
continues and g^, is updated every time that the particle en- 
ters element n, until the trajectory closes on itself for peri- 
odic Orbits or becomes dense in the tu-space. The only extra 
effort of this procedure compared to Schwarzschild's method 
is to perform the integral (|24p . In ijS.Sl we show that one can 
avoid this numerical Integration for separable modeis. 



3 GALERKIN WEIGHTING OF GOVERNING 
EQUATIONS 

This section implements a Bu bnov-Galerkin procedure 
l|Zienkiewicz. Tavlor fc Zhrj|2005l l to satisfy the governing 
equations of physical quantities (dependent variables) over 



individual elements in a weighted residual sense. As a re- 
sult, independent variables are eliminated from equations, 
leaving a System of algebraic equations between nodal val- 
ues of DF, density and velocity moments. The formulation 
is done in Cartesian coordinates and it should be modified 
for non-Cartesian ones (see iJS]for spherical Systems). 

3.1 Density and velocity moments 

Inside the nth element in the configuration space, equation 
((TJ reduces to 



Hn{x)pn{x) = / Hn{x)f{J) dv, 



(25) 



the presence of Hn {x) on the right side ensures that the Inte- 
gration is carried out over a phase subspace whose particles 
Visit the nth element and contribute to the density and ve- 
locity dispersion of that element. By substituting from (|16|l 
and p8|) into (|25|l we obtain 

M . 

H„{x)g„{x) ■ 6„ = J2H--(J) / Hn{x)f^{J) ■ P^ Av. (26) 

We now left-multiply this equation by da; g^{x) and inte- 
grale the result over the a;-domain to get 



(27) 



(28) 



G„ • 6„ = ^ / /da; Av Hm{J) 

TTl— 1 

X [_H'„(a;)gJ^(a;)-/„(J)] ■ p„ 



with 

Gn = 



H-a{x) ygl{x) ■ g„{x)^ da;. 



being an NdXNd constant matrix. The function H„{x) in the 
integrand of ([28} restricts the domain of Integration to the 
region occupied by the nth element. The integral in G„ can 
be done analytically should one use Interpolation functions 
of polynomial type. 

The transformation {x,v) — >■ {w,Xj is canonical and so 
the infinitesimal phase space volume da; dv can be replaced 
by dw dj. Using (|23|) with h — h ~ 0, equation (|27|l is 
transformed to 

ftn = f]^ /"/"d«;dJH„(J)e"='" 



xG„^ • [gfc(n,J)-/„(J)] p^ 



(29) 



It is obvious that only the term with k — (0, 0, 0) = con- 
tributes to the integral over the ly-space and equation 
reads 



&„ = ^ Fo(n,m) -p^, n = l,2,...,N, 



with 



(30) 



F4n,m) = {27vfJ dJ H^{J) G-' ■ [fo{n,J)-f^{J)]. (31) 

Repeating the above procedure for the functions Hn{x)u'{x) 
and Hn{x)T^-' (x) leads to 

M M 

4 = X] Uo(i,n,m) ■ p^, dt^ = ^ So(i, j,n,m) ■ p„, (32) 
© 0000 RAS, MNRAS 000, 000-000 
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where 



Uc = {2-Kf j dJ H^{J) G-^ ■ [gl{i,n,J)-f^{J)\ , (33) 
Se = i27TffdJH^{J)G-'-[f,ii,j,n,J)-f^iJ)]. (34) 



The constant element matrices Fe, Uc and Sc are of the 
size A^d X A/d and there are N x M of them. There are ad- 
ditional constraints associated with the element equations 
H3üp and (|32|) at a node shared by several elements, since a 
physical quantity must have the same value in the Galerkin 
projections of all those elements. In fact, one can intro- 
duce A't-dimensional column vectors b, Ci and dtj that con- 
tain all nodal densities and first- and second-order veloc- 
ity moments, and because nodes are shared the dimension 
Nt < N X Nd- Similarly, the nodal DFs constitute an Mt- 
dimensional column vector p where Mt < M x Md is the 
total number of distinct nodes in the J-space. Equation (|30p 
can thus be written as 



b=F p. 



(35) 



This matrix equation can be solved to yield the DF, as 
parametrized by its nodal values p. The rank of the ma- 
trix F is generally less than its dimension A't, and there are 
additional constraints that the DF should be non-negative, 
so in either Schwarzschild's method or the FEM the Solu- 
tion must be sought by linear or quadratic programming or 
some other optimization method (see i]3.4p . Once the Solu- 
tion is known, the nodal values of the Streaming velocity and 
velocity-dispersion tensor are obtained from equations p2p . 
which can be written as 



Ci = U(i)-p, dij = S{i, j) ■ p. 



(36) 



The A^t X Mt constant global matrices F, U and S are gen- 
erally dense. Assembling the element equations is a routine 
procedure in finite element analysis, and its logic is to use the 
continuity condition and eliminate repeated nodal quantities 
(like density and DF) from all element equations except one. 

3.2 Continuity and Jeans equations 

The accuracy of FEM modeis of equilibrium stellar Systems 
can be improved by adding additional constraints that en- 
sure that the continuity and Jeans equations pop and (llip 
are satisfied. Inside the nth element, these equations can be 
written 



H„(x)E|^-4 



0, 



(37) 



3 



dg-, 



9$ 



^"(=^) E P" • * = ^H4x)^g„ ■ 6„, i = 1, 2, 3. (38) 



Defining 



Hnix) 



9n ■ 

9$ 



dOn 



dxi 



dx, 



(39) 



Un(i) = - J H4x) — [gl ■ g,^] dx, (40) 

one obtains the Galerkin projections of p7|l and (138p as 

3 3 

J2'^i-ci = o, ^[u„(i)]-^-r„.d:f = &„. (41) 



We assemble these element equations to obtain the global 
forms 



^A(j).c,=0, ^B(j).d., = 6. 
Combining ((36]) and (|42ll leads to 



(42) 



J2 [ A(j) • U(j) ] ■ p = 0, ^ [ B(j) . S{i,j) ]-p^b. (43) 



j=i 



j=i 



We make some remarks. (i) The Solutions of the continuity 
and Jeans equations, whether the continuous versions (|lüp 
and pip or their FEM counterparts (|42p above, are gener- 
ally not unique. A notable exception is triaxial potentials 
of Stäckel form, in which the seco nd-order tensor t'^x) is 
diagonal in ellipsoidal coordinates (|van de Ven et al.ll2003l ). 
(ii) When using Co finite elements, as we do here, the mo- 
ments p, u^ , and r'-* (eqs.fTHU ^-^e continuous across element 
boundaries but generally their derivatives are not; however, 
since the right-hand sides of the continuity and Jeans equa- 
tions are continuous across boundaries, the combinations of 
derivatives of it' and r'-' that appear on the left-hand sides 
of these equations must also be continuous. (iii) Equations 
(|43p do not say that the mass and momentum flows into 
each element, through its boundaries, exactly balance their 
outflows (although the balance will become more and more 
accurate as th e number of no des increases). Only finite vol- 
ume methods (|LeVeauel2002l ') and conservative FEMs ensure 
the exact conservation of physical fluxes and this paper does 
not discuss those techniques. 



3.3 Separable modeis 

The computation of the element matrices Fe, Uo and So is 
accelerated when the Hamilton-Jacobi equation is separable 
for the Potential $(a;). In such a circumstance, the velocity 
vector can be expressed as v{x,J), which implies 



dv = Q{x, J) dJ, Q{x, J) 



d{vi,V2,V3) 



(44) 



9(Ji, J2, J3)' 

This allows us to bypass the costly Integration of orbit 
equations needed for calculating the Fourier coefRcients gg. 
Defining the matrix 

P{x, J) = H^{x)H^{J) Q{x, J) G-i ■ [glix) ■ /„(J)] , (45) 

the element matrices are computed from 



f I P{x, J) dx dJ, 



Vi{x,J) P{x,J) dx d J, 



/ / Vi{x,J) Vj{x,J) P{x,J) da; d J. 



(46) 
(47) 



(48) 



The Integrals over the x and J subdomains can be evaluated 
using Gaussian quadratures. 

In separable modeis the turning points of orbits and 
their shapes in the configuration space are known. Therefore, 
the null Integrals in the element matrices can be avoided, and 
the computational effort is reduced, by a priori Identification 
of the J-subspaces whose orbits never enter a selected ele- 
ment in the a;-space. In fact, the Information related to the 
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passage of an orbit through a given element is carried by the 
function Hm{J)Hn{x), and before evaluating the Integrals 
we can exclude all {m,n) pairs that give H„{x)Hm{J) = 0. 
For separable modeis in non-Cartesian coordinates, the ve- 
locity components in the Jacobian Q are replaced by gener- 
alized momenta, and the matrix P is divided by the product 
of metric coefficients. 



3.4 Linear and quadratic programming 

The size of the unknown vector p is not necessarily, or usu- 
ally, equal to the total number of constraints. Even if it 
were, the Solution vector would not necessarily fulfiU the 
requirement that the DF must be non-negative. We there- 
fore employ either linear programm ing (LP) or quadratic 
Programming (QP; iGill et al.lll98ll ) and search for pi, the 
components of p, by minimizing the objective function 



Mt , A/t Aft 



pipi' 



(49) 



;=i i' 



with Will ~ for LP. The minimization is subject to the 
inequalities p; > (for 1=1, 2, . . . , Mt) and the equality con- 
straints p5|) . If we also demand satisfaction of the continuity 
and Jeans equation these are supplemented by the equality 
constraints (|43p . The QP routines begin from a vector Pq 
that satisfies the imposed constraints with a tolerance of e/, 
then proceed to minimize J7. The vector Pq is usually called 
a feasible Solution and e/ is the feasibility tolerance; the lat- 
ter must be greater than the computational accuracy of the 
variables involved in the constraints. 

The weights Wui are chosen based on the desired at- 
tributes of the model, such as blas toward radial or tangen- 
tial Orbits, maximization of a quadratic entropy, or fitting 
to specified data. For example, if a series of observables Oa 
are linear functions of the DF, 



A/t 



Oa=2_,OalPl, a = 1, 



(50) 



and they are observed to have values Oq with observational 
errors aa, then a suitable objective function is specified by 

G = -^ö„0.,, Wa' = J2^^^- (51) 

— — ö"^ 

The LP and QP algorithms we have used can stall at weak 
local minima or "dead points" . Whenever this happens, we 
perturb the Solution and restart the algorithm. 



4 DERIVATION OF SCHWARZSCHILD'S 
METHOD 

It is now straightforward to show that Schwarzschild's or- 
bit superposition method is a subclass of FEM. We assume 
for simplicity that the potential is integrable so the orbits 
are regulär. The orbit library in Schwarzschild's method is 
coUected by sampling over the space of initial conditions. 
For regulär orbits there is a one to one and invertible map 
between {xo, vo) and [w{0), J], and Schwarzschild's DFs will 



have the foUowing form (|Vandervoortlll984l ) 



M 



^^"^ " 72^ ^ P™.'5(J- Jm), 



(52) 



where &{■ ■ ■) is the Dirac delta function and pm is the dis- 
crete DF associated with an orbit of the action vector Jm- 
Equation (|52|l is derived from (118p by shrinking the elements 
in the action space to zero size. The total mass of the galaxy 
is computed from 



M 



f{J) dJdw, 



which is combined with (I52II to obtain the constraint 



E P™ ^ 1- 



(53) 



(54) 



Schwarzschild assumes a uniform density p„ inside the 
nth element in configuration space. This implies that there is 
one node per element (A^d = 1) and that the vector function 
g„{x) reduces to a scalar constant, g„ = 1. Equation (|12|l 
then reduces to 



P(^) = ^H„{x) 



(55) 



The matrix Gn is now a Single number Vn, which is the 
volume of the nth element. The quantity A4n=Gn ■ bn=VnPn 
will thus be the mass inside the nth element. We Substitute 
((52)l and ([SSj into (|3ÜJ and obtain 



Mn = AI ^ Pm go{n, Jm), 



with the zeroth-order Fourier coefficient given by 
go{n, J,r ' 



(27r)3 



Hn[x{w,Jm)] AW. 



(56) 



(57) 



According to time averages theorem (|Binnev fc Tremaind 
2008), the quantity on the right hand side of (157p is the 
fraction of time tn{Jm) that an orbit of action Jm spends in 
the nth spatial element. Consequently, we obtain 



X„ = X E *"('^" 



Pm, 



(58) 



which is Schwarzschild's equation. 

The approach described in tj3.2l to incorporate con- 
straints based on the continuity and Jeans equations into 
FEM modeis does not work for Schwarzschild's method: be- 
cause the interpolating functions g are constants, the ma- 
trices T^ defined in equation (|39p are zero so the first of 
equations (|4ip is trivially satisfied and the second has no 
Solution. Physically, the stress tensor t'Ü' is constant within 
an element so there is no divergence in the momentum flux to 
balance the gravitational force per unit volume —pd^/dxi. 

The failure of Schwarzschild's method to satisfy the 
Jeans equations within an element does not imply that the 
method is incorrect. Indeed, the method must satisfy the 
Jeans equations on larger scales because the assumed dis- 
crete DF (|52p depends only on the actions and hence must 
satisfy the CBE, and the Jeans equations are moments of the 
CBE. The correct Statement is that Schwarzschild's method 
satisfies the Jeans equations approximately if we calculate 
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gradients of the stress tensor between adjacent elements and 
match their sum to — pV"!" at the center of an dement. In 
this process, which we carry out for spherical Systems in JSl 
one must appropriately handle partial derivatives because el- 
ements in the configuration space are not usually separated 
uniformly. 



5 SPHERICALLY SYMMETRIC SYSTEMS 

In spherical Systems one can use simple shell elements (ring 
elements for axisymmetric disks). Moreover, the distance of 
particles from the centre is represented as the Fourier se- 
ries of the radial angle wr only. This remarkably simplifies 
the calculation of the vectorial function Qq should one de- 
cide to compute the dement matrices from (|31|) . (|33p and 
((34]). Consider the spherical polar coordinates (r,6,(j>) and 
the corresponding velocities {vr,ve,v^) where r is the radial 
distance from the centre, 9 is the co-latitude and 4> is the 
azimuthal angle. The density of a spherical System is a func- 
tion of r, its first-order velocity moment {vr) is zero, and the 
foUowing relations hold for second-order velocity moments: 



{vi) = 2{v'l) = 2{vl), (VrV^) = (VrVe) = 0, 



(59) 



where Vt = Vg +v'^ = 1? jr^ with L being the magnitude of 
angular momentum vector. We confine ourselves to modeis 
with {v^} — {v0) — {v^ve)—O. The continuity equation (mass 
conservation) is trivially satisfied for such a System. The 
elements of the stress tensor are 

r"' = p{v1), r" = p{v1), r®* = p{vl), t''"*' = p{vl), (60) 

and the Jeans equation in the radial direction reads 



^iZ. + i ('2r'''' - r") = _ — 
Ar r ^ ' dr 



(61) 



The other two equations, in the 9- and (^-direction, do not 
provide further Information. 

We consider a mesh of A'^ concentric shell elements and 
dehne the nth dement by its inner and outer radii r„ and 
r„+i, respectively. We then approximate the density and 
velocity moments as 



n = l 

JV 

r'-'-(r) = ^i/„(r)gjr).<:, 

JV 

r"(r) = Y.HÄr)grXr)-€. 



(62) 
(63) 
(64) 



The continuity conditions at the boundaries of adjacent ele- 
ments are &i,{n-i-i) = ^iVd,«, rfiy{„+i) = <^iVd,n and d"(„^i) = 
d% n- Let US now Substitute from equations (|62|| - (|64p into 
(|6H) and derive its Galerkin protection as 



^ V-i . [T; • S7(n, m) - n ■ S**(n, m)] ■ p„ = 6„, (65) 



for n = 1, 2, . . . , TV. The matrices Se"^ and Se' are determined 
from (|48|) by replacing ViVj with w^ and u?, respectively, and 



we have 

' n 

T' 



H„{r) 
= - f H, 



^" dr 
1 T 

d$ r T 

\9n 



r'dr + 2Tl,, 



(") d^l»"-^"] 



r dr. 



(66) 
(67) 
(68) 



All other equations of ij3]can be directly applied to spherical 
Systems using the foUowing substitutions: 

, 2 , , ^, ^,^ AtvL dEdL 

(69) 



(70) 



If we apply the FEM without Jeans equation constraints 
we must satisfy the linear constraint equations b — F ■ p (eq. 
I35|l . If we include the Jeans equation constraints we assemble 
equations (I65|l to a global form T ■ p — b and combine this 
with 6 = F • p to give 





r^Vr 


whf 


3re E is the orbital energy of particles: 


E = 


1 2 i' ^, , 
= 2"' + 2r2 +*('■)• 



F) • p = 0. 



(71) 



In the Co finite element formulation, it is difhcult to 
construct a function (here the stress components) and its 
derivatives with the same accuracy. Therefore, we replace 
the equality constraints (|71|l with the weaker inequality con- 
straints 



< e„ < En 



(72) 



where e„ are the normalised components of the residual vec- 
tor (T — F) ■ p, defined as 



, A/t 

= T-^(rn!-Ki)P^ n=l,2,...,iVt, 



and minimize J by setting 



JVt 



^^-E 



T 



bn 



(73) 



(74) 



Here T^i and F^i are the elements of T and F, respectively. 
The value of fmax cannot be arbitrarily small: at large radii 
the magnitudes of the stresses become comparable to the 
numerical errors, and the Jeans equations are correspond- 
ingly less accurate. In our calculations, we have been able 
to secure convergence with tmax as small as ^/ej over a wide 
radial ränge. 

It is worth deriving the weighted residual form of the 
Jeans equation for discrete Schwarzschild modeis, to inves- 
tigate whether applying this as a constraint improves the 
accuracy of these modeis. For g„ = 1, the density and sec- 
ond velocity moments are constant inside each element and 
the vectors b„, d^ and d^* in (|62|| - (|64p are replaced by the 
scalars p„, r,^"^, and r", respectively. We rewrite (|6ip as 

r^ dr r ar 

We obtain the Galerkin projection through multiplying the 
differential equation (|75p by H„{r)r dr and integrating over 
the nth spatial element. For the first term this procedure 
gives F{rn+i) — F{r„) where F{r) — r'^T^''{r). Since r'"'" is 
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discontinuous at the element boundaries we must make some 
arbitrary choice; after some experimentation we have found 
that the best convergence is obtained by taking r'^^irn) to be 
T,i!li, that is, the value of the stress in the element interior 
to the boundary. Then the discretized Jeans equation is 



+ 



-T„ = Pr. 



'" + ld"I> 2 

— — r Ar. (76) 
dr 



This difference equation is not necessarily satisfied by an 
optimization procedure that attempts to fit observables us- 
ing a DF of the form (|52|l . We note that r^li and r^"^ are 
normal stresses exerted on the nth shell element at its inner 
and outer boundaries. For sufficiently thin elements when 
(Ar„)^ < r„Ar„, one may replace r^+i by r^ + 2r„Ar„ 



and write (|76|l as 



Ar„ 



+ 



2r: 



d£ 

dr 



(77) 



which is the discrete counterpart of (|61l) obtained using a 
backward difference scheme. 



6 EXAMPLES 

We illustrate the Performance of FEM for spherically Sym- 
metrie S ystems by c onstru cting ergodic and anisotropic DFs 
for the iHernguistl (| 19901 ) model. Each orbit is character- 
ized by its majcimum and minimum galactocentric distances 
rmax(J) and rmin(J). We define 



a{J) 



+ rn 



e(J) 



: + rn 



(78) 



and the finite-element mesh is generated in (a, e)-space in- 
stead of action space. 

The potential-density pair for the Hernquist model is 
given by 



$(r) 



1 



1 + r' 



p{r) 



27rr(l + r-)3 



(79) 



Since p diverges toward the centre, and to resolve the be- 
havior of functions in the central regions, the nodes of our A'' 
shell elements are distributed with equal logarithmic spac- 
ing, using the power law 



= 10 



--1l+aiy(n,N) 



y{n,N) = 



+ 



n-1 



2(iV + l) iV+1 



n — 1, 



.,N + 1. 



(80) 



We use simple double-node elements with A^d — 2 (no inte- 
rior nodes) and linear interpolating functions in the radial 
direction, 



9jr)=[^{l~f) i(l + f)] 



2(r - r„) 



1. (81) 



Having the grid Information and Interpolation functions, the 
matrices G„, T^, T^, and V„ can be calculated. 

The mesh in the two-dimensional (a, e)-space consists 
of M = Ma X Me rectangular elements, each with Md = 4 
nodes (Figure[2}. For the mth element sitting in the jth row 
and ith column of the mesh, the local coordinates (^, rj) are 



e-i — 



7+1 



Figure 2. Two dimensional finite element mesh in the (a, e)- 
spacc. Local variables (§,»7) vary in the interval [— l,-)-l] and the 
centre of element is located at (^, rj) = (0, 0). 



defined as 

2(a — fli 



^ 



V = 



0,1+1 — a.i 
2(e-e,) 



ej+i 



1, i = i,2,. 



,M„ (82) 

,Me, (83) 



where m, = (i — 1) x Me -f j and 



10 



-72+a2!/(«'.Ma) 



e,, =y(/,Me). (84) 

The minimum eccentricity in our grid is ei = i/(A/e + 1); 
we avoid zero-eccentricity orbits because of the singularity of 
the Jacobian Qix, J) at the circular orbit boundary of the 
action space. We do not have exactly radial orbits in our 
modeis as the maximum eccentricity in our grid is eM^+i = 
1 - i/(Me + 1). 

The four nodal coordinates of each element are given 



by 

(C3,r;3) = (+1,-1) 



(6,^2) = (-l,+l), 

(e4,»74) = ( + l,+l). 



(85) 



The DF at (^^,77^) is indexed by k and the foUowing Inter- 
polation functions are used inside the mth element 



fkM^,v) = ]ii + ^kO{i + rikV), k = 1,2,3,4. 



(86) 



These are smooth quadratic functions that behave linearly 
along element boundaries. DFs at the common nodes of ad- 
jacent elements must be equal in order to build a continu- 
ous /(J). This condition is taken into account in assembling 
the global matrices F, U and S. Our experiments show that 
FEM is not highly sensitive to the parameters 7^ and a^ 
(i — 1, 2) but they should be chosen so that at least one or- 
bit passes through each element. The cost of computations 
is remarkably reduced by minimizing the size M of the grid 
in action space while keeping the number of constraints con- 
stant. In such conditions, securing a feasible Solution p^ be- 
comes harder though the choices 71 = 72 and ai = 02 are 
often helpful when the same element nodes are used in the 
r and a Spaces. The reason is that a Solution p — F~^ ■ b 
always exists in the limit of a DF composed of circular or- 
bits, / — /o(a)(5(e^), and one can imagine smooth DFs of 
the form 

/=(l-Q),fo(a)5(e') + Q/i(a,e), < a < 1, (87) 

which are found by the optimizer through varying a and /i . 
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6.1 Ergodic distribution functions 

We begin our numerical experiments by constructing ergodic 
DFs; these give an isotropic stress tensor with 



2r'' 



f2S'' 



S") ■ p = 0. 



Here S^"^ and S" are Nt x Mt matrices assembled frorn S" 
and SJ*", respectively. We solve this particular problem by 
linear programming (LP). Our first FEM model has TV = 50 
Shell elements in the configuration space and a mesh of 
Ma X Me = 40 X 40 elements in the (a, e)-space. The grid 
points have been obtained by setting ai = 0:2 = 3 and 
71 = 72 = 2. The innermost grid point in the configuration 
Space is at r\ = 0.0107 and the outermost at rjv+i = 9.345 
where the Hernquist model density is 1.54 x 10~^. The inner- 
most and outermost grid points in the a-direction are located 
at a\ — 0.0109 and om^+i = 9.1921. The total number of 
unknown nodal DFs is Mt = 1681. For A^ linear elements 
with A'd = 2 nodes per element, we have A'^x (A^d — l)-l-l ~ 51 
equality constraints to build p from b = f ■ p, and 51 equal- 
ity constraints to impose the isotropy condition (|88p . There 
are also 51 inequalities of the type (|72p when the Jeans con- 
straints are present. The weights of all nodal DFs in the 
objective function J = "^^ Cipi are assumed to be equal: 
C; = 1 (Z = 1, 2, . . . , Mt) in the absence of Jeans constraints, 
which corresponds to minimizing the sum of the values of the 
DF at all the nodes. Experiments with other choices for the 
Ci suggest that our Solutions are not sensitive to this choice, 
which is to be expected since the ergodic DF for a spherical 
System with a given density and potential is unique. 

We also construct a model with the parameters A'' = 50, 
Ol = «2 = 3 and 71 = 72 = 2 as in the first model, but with 
a coarser grid of Ma x Me = 30 x 30 in the (a, e)-space. 
These give ai = 0.0112, öm^+i = 8.9457 and Mt = 961. 

In all of our runs with and without Jeans equation con- 
straints, p{r) is successfuUy reconstructed with a fractional 
error < 10~*, and the anisotropy parameter ß{r) = 1 — 
ir^/r'"'" is zero to within the feasibility tolerance e/ = 10~*. 
The fractional accuracy in satisfying the Jeans equations 
is controlled by the parameter tmax (eg. 172p . We initialize 
emax to 0(10"*) for r < 3.2 and 0(10"^) at larger radü 
where the magnitude of r'"'" becomes comparable with the 
discretization errors, which are greater than the feasibility 
tolerance by several Orders of magnitude. At some nodes the 
prescribed emax may be too small to allow for a reasonable 
optimal Solution. We tolerate constraint violations of up to 
5% at individual nodes should the RMS of tn (over all nodes) 
beof 0(10"^). 

Figure[3l3 displays the computed density p{r), radial ve- 
locity dispersion ö-,(r), and fractional error e^ = 1 — ar/a^, 
where the exact dispersion a'^. comes from equation (10) in 
iHernauisd (Il990l ). The FEM Solution constrained by equa- 
tion (|88p but not the Jeans equation constraint (|7ip exhibits 
an error of 3% in the outermost element, and a systematic 
drift from a-^{r) in the central regions, amounting to a 1.5% 
error for r < 0.03. If we had not any Information about 
the exact dispersion profile, the computed ar was smooth 
enough to be accepted as a Solution. When we add the Jeans 
equation constraint the mean error is reduced by a factor of 
2.5, typically to < 0.3%. The errors of up to 0.5% near the in- 
ner and outer boundaries are due to FEM discretization and 
model truncation, and hence are not improved by adding the 



Jeans equation constraint; these can be suppressed by using 
boundary elements to cover the currently neglected ranges 
[0, ri) and (riv+1,00). Adding more spatial elements is not 
helpful beyond the radius at which the stresses become as 
small as the discretization errors. 

To compare FEM with Schwarzschild's method, we 
build Schwarzschild modeis using the same TV = 50 shells 
in configuration space, with a discrete DF (eq. 152p that is 
non-zero only at actions Jm given by the nodes of the mesh 
defined in (|84p . The equality constraints in the LP routine 
consist of Schwarzschild's equation (|58p and a variant of 
the isotropy constraints H88p . We use the same tolerances 
as in the FEM modeis. The profiles of p, ar and €„ in our 
Schwarzschild modeis are shown in Figure[3f) for two grids, 
Ma X Me = 30 X 30 and 40 x 40. It is seen that the dispersion 
error e^ is as large as 10%, about five times larger than in the 
FEM method; the radial fluctuations in e^ are also larger. 
We remark that the model with the smaller orbit library 
or action-space grid (Afa x Afe = 30 x 30) is more accu- 
rate, which highlights the fact that Schwarzschild's method 
is sensitive to the locations Jm of the delta functions in 
the action space. Special orbit sampling strategies can re- 
duce this sensitivity and enhance the accuracy of computa- 
tions (jThomas et al. 20041). Note that this sensitivity is not 
present in FEM; Figure[3li shows FEM modeis with 30 x 30 
and 40 x 40 grids and the error is generally smaller with the 
larger grid, as expected. 

Adding the constraints (|76p to the optimization pro- 
cedure brings a dramatic change for Schwarzschild modeis. 
For the model with Ma x Me = 40 x 40, the curve of ar 
and its fractional error are astonishingly smooth when the 
Jeans constraints are imposed, although the rms error is not 
substantially reduced. By following this procedure, the sen- 
sitivity to the sampling of orbits (or choosing the location 
of delta functions) disappears. Our experiments with modeis 
constrained by the Jeans equation show that the deviation 
between the exact and computed curves of ar is of 0{Ar„). 
An exact match is anticipated in the limit of Ar„ — >■ 0, but 
convergence to the correct Solution will be slow and costly. 
In contrast, FEM converges at a rate 0{h^'^^) where h is the 
element size (in either action space or configuration space), 
and p is the order of the interpolating polynomial (p = 1 in 
our case). Thus FEM should, and does, converge faster than 
Schwarzschild's method by one order in the element size. 



6.2 Anisotropie distribution functions 

In this subsection we use FEM and QP optimization to con- 
struct anisotropic DFs. We begin with a model constructed 
without Jeans equation constraints (e„ varies freely) having 
TV = 70 spatial sheU elements and a mesh Ma x Me = 70 x 70 
in (a, e)-space. We set 71 = 72 = —1.5 and «1 = 02 = 2.5, 
which correspond to ri = ai = 0.0329 and rjv+i = aM„ + i ~ 
9.6027. The weights in the objective function (|49p are chosen 
as Ci = and 

Wa'{J) ^ Sa>[l-e\ji)], l ^ (i' - 1){M, + 1) + f , (89) 

i' = l,2,...,Ma + h 
/ = l,2,...,Me-H, 

which is designed to minimize the population of low- 
eccentricity orbits. Note that Wui is positive definite so there 
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Figure 3. The computed profiles of the density p, radial velocity dispersion (Tr, and fractional dispersion error f.a for the Hernquist model 
with ergodic DF. In all FEM and Schwarzschild experiments, the anisotropy parameter /3(r) is zero over elements with an accuracy of 
10~*. The exact values of p and a^ are shown by filled circles. (a) FEM results for N = 50 Shells in the configuration space. Three modeis 
are shown, although they are indistinguishable in the top panel and almost indistinguishable in the iniddle panel: Ma x Me = 40 X 40 
elements in the (a, e)-space, no Jeans equation constraints (dashed line); Ma x Me = 40 X 40 with Jeans constraints (solid line), and 
Ma X Me = 30 X 30 with Jeans constraints (dot-dashed line). (b) Schwarzschild modeis with N = 50 Shells. The three modeis are: 
Ma X Me = 30 X 30 elements in (a, e)-space, no Jeans equation constraints (dashed line); Ma x Me = 40 X 40, no Jeans constraints 
(dot-dashed line); Ma x Me = 40 X 40 with Jeans constraints (solid line). Note that the vertical scales in the two bottom panels are 
different. 



is a Single global minimum of the objective function. The [S] For this model, the RMS error 

node / in the (a, e)-space is located at {aii,eji) where a^/ 



and Cji are computed from (184p . We thus have e(J;) — e^. . 
There are {Ma + 1) x {Me + 1) = 5041 unknown compo- 
nents of p, and 71 equality constraints that correspond to 
b = f ■ p. The QP routine converges and finds the global 
minimum J = 3.218 for the objective function. Figure |4]a 
shows the isocontours of the computed DF. The smoothness 
of the DF is evident. The narrow curved feature is inher- 
ited from a feasible Solution: the QP method smooths the 
distribution around a feasible Solution pg that satisfies Prob- 
lem constraints. Since the number of non-zero components 
of Po is much less than Mt, the subdomain covered by that 
feasible Solution in the action space shows up as a distinct 
feature. 



To probe whether the Jeans equation is satisfied in this 
model, we have calculated the distribution of the errors e« 
in the Jeans equation (|73|) . and we display these in Figure 



1 ^' 



A^t 



1/2 



(90) 



is e = 0.226. It is evident that the differential Jeans con- 
straints (|65p have been violated by a large margin in the 
innermost element, and this problem bleeds over into the 
Stresses in nearby elements, out to r ~ 0.3. 

We now include the differential Jeans equation con- 
straint (|71|l in the Solution procedure. Since the number of 
constraints has been increased, the QP algorithm yields a 
larger objective, J = 54.607. We find e = 0.056 which is 
four times smaller than the RMS error in the unconstrained 
model. The relatively large errors in the outermost elements 
are due to discretization errors, just as in the case of the 
ergodic Solutions of the preceding subsection, and could be 
corrected by a proper boundary element. 

The computed DF (displayed in Figure Hf) is now less 
smooth and has developed several narrow eccentricity spikes, 
the strongest of which is at e ~ 0.75. The centroid of the 
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Figure 4. Isocontours of logig{f) for radially biased modeis unconstrained (panel a) and constrained (panel h) by the differential Jeans 
equation II6II I. Both modcls have 70 X 70 elements in the (a, e)-space and 70 shell elements in the configuration space. 
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Figure 5. The residuals e„ and radial velocity dispersion ar = 
(i)^)^'^ for modeis constrained (solid lines) and unconstrained 
(dotted lines) by the differential Jeans equation JüTJ. The node 
number n and its corresponding radial distance r„ are related 
through log]^Q(rn) = — 1.5+2. 5y(n, TV). The error in the innermost 
bin is off the scale of the lower panel, ei = 1.828. The larger 
errors in the outermost elements arise because the magnitude of 
the stress r'''" is comparable to the discretization errors for r > 8. 



DF has also shifted to smaller a. The narrow curved feature 
has also shifted, because the number of constraints has been 
doubled and a new feasible Solution has emerged. There are 
several isolated small rectangular regions with non-zero /; 
these would have been delta functions in Schwarzschild's 



approach, but now occupy subdomains containing at least 
four elements. The continuity and difFerentiability of / is 
clear even in such isolated subdomains. The radial velocity 
dispersions of the two Solutions (with and without Jeans 
equation constraints) are compared in Figure [S] 

Since some wiggles exist in the dispersion ar of the 
constrained model and the magnitude of J is substantially 
larger than in the unconstrained model, we suspect that QP 
has not been able to find a global minimum correspond- 
ing to a reasonably smooth DF. This may have occurred 
because our strict e„ — > 0"'' constraints have shielded the 
global minimum. To locate the global minimum, one can 
identify unnecessarily small values of e„ and ease the cor- 
responding constraints that may have strongly constrained 
dr'^''/dr in a p-subspace where this gradient is actually in- 
accurate due to errors in the computed t'''". An alternative 
way of employing Jeans constraints, which gives smoother 
DFs, is discussed below. 

We integrale equation (I61|l to obtain 



rr ar 



r p——dr, 
ar 



(91) 



in which we have imposed the boundary condition that the 
Stresses vanish at infinity. Substituting from (|62p - (|64p into 
(|9ip and performing Integrals over individual elements, gives 



JV 



N 



J2 Hn{ryg„ -dZ^Yl H4r) (3" ■ d" + s! ' bn) , (92) 

n— 1 n^l 

where 

,d<l> 



Sn ('■■) = / rgjr) dr, gl{r) 



Ar 



gjr) dr. 



Multiplying H92p by drH„{r)g'^{r) and integrating over the 
r-domain leaves us with 

M 

^ V„ ■ I^Gn- 57(71,771) -t„ ■ S"(n,m)J -p^ = 6„, (93) 

TTi — 1 
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Figure 6. Top: Isocontours of log]^Q(/) for radially biased modeis 
constrained by the integral Jeans equation. There are 70 X 70 ele- 
ments in the (a, e)-space and 70 shell elements in the configuration 
Space. Bottom: radial velocity dispersion ar and the normaliscd 
error e„. The solid and dashed lines correspond to models with 
and without integral Jeans equation constraints H93| l , respectively. 
The DF of the unconstrained niodel is displayed in Fig. |4]2. 



ior n= 1,2, 



, N with 



\l^ = j H,,{r)gl-gtdr, fi = J H„{r)gl ■ g'^ dr. (94) 



One can assemble equations (|93p into a global form T p = b 
and foUow the procedure of !j5]by replacing T with T in equa- 
tions (|7ip - (|74p . The advantage of these integral Jeans con- 
straints (|93|l over their differential counterparts (|65|l is that 
they do not involve derivatives of the interpolating functions 
and r"". 

We have assumed the same niodel properties as in Fig- 
ure |4] and constructed a smooth DF (top panel in Fig. [6]) 



using the integral Jeans equation. To impose the constraints 
T ■ p = 6 we have minimized the new residuals 



1 "' 



PI, n- 



1,2, 



,iVt. 



(95) 



using only the lower zero bound on them: e„ > 0. Figure [6] 
displays ar{r) and the Variation of £„ for models that are un- 
constrained and constrained by the integral Jeans equation. 
The close agreement between the radial dispersion curves of 
the unconstrained and constrained models shows that the 
DF of the unconstrained model, which is identical to the 
model of Figure \^, is accurate and smooth enough to sat- 
isfy the integral Jeans equation. For the constrained model, 
we have found J = 3.453 which is close to the unconstrained 
minimum J = 3.218, indicating convergence to the global 
minimum. Smaller values of e„ in the central regions of the 
unconstrained model, compared to what we displayed in Fig- 
ure[5]for e„, show that the DF constructed by FEM satisfies 
the integral Jeans constraints more accurately than differen- 
tial ones. The correction e„ — ^ 0^ in the constrained model 
yields satisfactory results in the central regions, but again, 
the outermost dement shows a large error because of the 
model truncation. 

We are also able to construct anisotropic models with 
prescribed anisotropy profiles. In terms of /3(r), the Jeans 
equation reads 



dr"- 



2ß{r) 



d$ 



, , . , , . (96) 

dr r dr 

We seek a model with ß{r) = —1/(1 4- r) that has a tangen- 

tially biased core and becomes Isotropie as r — ^ oo. For the 

Hernquist model, this yields the exact Solution of equation 

dMll: 



6 In 



27r(l-^r)2 

l+r\ {l + 2r){6r^ + 6r - 



1) 



2r2(H-r)2 J- '^^''^ 

We now foUow the procedure of SJS] and project the equation 
2(1 — ß)T^^ = r" on the p-space. One can verify that 



E 



IT 



S7'{n,Tn) - -S"(n,m) 



Pm = 0, 



(98) 



holds for n = 1, 2, . . . , A^ where I is the identity matrix of 
dimension A'^d x A'd and 



ff„(r)/3(r)G; 



T 

g g. 



n] t'' 



dr. 



(99) 



After assembling the System of new constraints (|98|) into 
a global form, we ran our FEM code and compared its out- 
come for Gr with that of equation (|97|l . Figure [7] illustrates 
our results for A'^ = 50 spatial shells and Ma x Me = 30 x 30 
elements in the (a, e)-space. As in the ergodic case, there 
are distinguishable differences between models constrained 
and unconstrained by the Jeans equation. We did our cal- 
culations using both LP and QP, and could recover theo- 
retical curves of p and ß with the fractional accuracy 10~* 
in both approaches. The QP results were not sensitive to 
the choice of Wui, but we worked with Wui = 5iiie{Ji) 
{l = 1,2, . . . , Mt), which is more consistent with a tangential 
core. For the LP without Jeans constraints, choosing Ci = 1 
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Figure 7. Density, anisotropy parameter, radial velocity disper- 
sion and fractional error Eo- for a model with /3(r) = —1/(1 + r). 
The exact values of quantities have been shown by filled circles. 
The solid and dashed lines correspond to models with and without 
differential Jeans equation constraints II65I I. respectively. 



always gave accurate results. The results of LP and QP were 
almost identical. 



prove the convergence of FEM-constructed equilibrium DFs, 
our numerical experiments with the Hernquist model show 
convergence to the exact values of observables as the el- 
ement sizes are decreased in a uniform logarithmic mesh. 
Nevertheless, a major challenge in the convergence analysis 
of equilibrium models constructed by FEM is to understand 
the role of constrained optimization routines and why they 
appear to give non-unique models in some cases such as Fig- 
ure [S] 

The disadvantage of FEM codes is that they are more 
time-consuming to write and to run than Schwarzschild 
Codes — but not by a large factor. The main difference is 
in the subroutine where the element matrices are computed 
and assembled, and in our implementation this subroutine 
contains fewer than twice as many Statements. The time re- 
quired for optimization is the same, since the matrices have 
the same size in both methods, but the construction of the 
matrices takes longer in FEM. Overall, the calculations used 
to produce an FEM model usually take 3-5 times as long 
to run as the calculations for the analogous Schwarzschild 
model. 

There are at least two reasons why the higher accuracy 
provided by FEM is likely to become more important in the 
future. The first is that the quality of kinematic and Pho- 
tometrie observations of galaxies is growing, in particular 
because of integral-field spectrographs on large telescopes, 
and higher quality data demand more accurate dynamical 
models. The second is that the State of the art has gradu- 
ally progressed from spherical models to axisymmetric and 
triajcial ones; the higher dimensionality of triaxial models 
demands much larger grids in configuration space and this 
in turn calls for numerical methods that converge as a higher 
power of the characteristic scale of the grid Clements. 

It is also possible to apply FEM to time-dependent stel- 
lar Systems, such as barred and spiral galaxies. Since the 
evaluation of the element integrals does not depend on the 
nodal variables, the extra computational cost of unsteady 
Problems is associated only with the integration of the evo- 
lutionary equations in the time domain. Ijalalil (|2010|) has 
used FEM to study the linear stability of razor-thin disk 
galaxies and the response of such galaxies to external per- 
turbations such as satellite galaxies. 



7 DISCUSSION 

We have demonstrated that Schwarzschild's method for con- 
structing stellar Systems can be regarded as a special case of 
finite element methods (FEM), and that FEM can be used 
to construct stellar Systems that are more accurate approx- 
imations to the coUisionless Boltzmann equation for given 
grids in configuration and action space. We have also shown 
that the accuracy of both Schwarzschild and FEM models 
can be substantially improved by incorporating the Jeans 
equations as explicit additional constraints. 

There are rigorous mathematical methods for prov- 
ing the convergence of Co FEM schemes to continuous 
and smooth Solutions of initial and boundary value Prob- 
lems described by ordinary and partial differential equa- 
tions l|Szab6 fc Babuskal ll991h . Such analyscs have also been 
carrie d out for eigenvalue problems ( B abuska fc Rheinboldtl 
1 1978h and adaptiv e FEMs in two dimensional problems 
( Morin et al.l 12002 ). Although we did not mathematically 
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